Build a phyloseq object

Let’s start by reading your data:

# Import libraries:
# Load libraries
# Note: commented-out libraries are experimental libraries for other functions not used in this final analysis file
library(limma)
#library(DESeq2)
library(edgeR)
#library(Rsubread)
library(sns)

Package: sns, Version: 1.1.2
Metropolis-Hastings MCMC using Stochastic Newton Sampler
Scientific Computing Group, Sentrana Inc. &
Imperial College London

#library(Glimma)
library(data.table)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.14.2 using 20 threads (see ?getDTthreads).  Latest news: r-datatable.com

library(tidyr)
#library(magrittr)
library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6     ✔ dplyr   1.0.9
✔ tibble  3.1.7     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.2
✔ purrr   0.3.4     
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::between()   masks data.table::between()
✖ dplyr::filter()    masks stats::filter()
✖ dplyr::first()     masks data.table::first()
✖ dplyr::lag()       masks stats::lag()
✖ dplyr::last()      masks data.table::last()
✖ purrr::transpose() masks data.table::transpose()

#library(sjPlot)
#library(sjlabelled)
library(sjmisc)

Learn more about sjmisc with 'browseVignettes("sjmisc")'.

Attaching package: ‘sjmisc’

The following object is masked from ‘package:purrr’:

    is_empty

The following object is masked from ‘package:tibble’:

    add_case

The following object is masked from ‘package:tidyr’:

    replace_na

library(ggplot2)
#library(effects)
#library(broom)
#library(GGally)
#library(emmeans)
#library(ggpubr)
#theme_set(theme_pubr())
library(corncob)
library(phyloseq)

# Import data:
otus <- read.csv("gut_full.csv", check.names=F)
taxa <- read.csv("taxa.csv",sep='\t')
otus$sex <- factor(otus$sex) #factorize sex
otus$vendor_dashboard <- factor(otus$vendor_dashboard)
otus$bowel <- factor(as.numeric(factor(otus$bowel, levels = c(1,2,3,4), labels = c(1,2,3,4))))
otus <- otus[ which(otus$vendor_dashboard == "Second Genome" | otus$vendor_dashboard == "research-microbiome"),] # keep only data where vendor is explicit
otus$vendor_dashboard = str_replace_all(otus$vendor_dashboard,"research-microbiome","DNA Genotek")
otus$vendor_dashboard <- factor(otus$vendor_dashboard) # factorize vendor
#otus[ which(otus$vendor_dashboard == "Second Genome"),]
otus <- within(otus, bowel <- relevel(bowel, ref = 3))

# Pre-processing and checking validity of data
# Algorithms provided by Christian Diener, PhD:
dat <- otus
otus <- otus[!duplicated(otus$public_client_id), ]
genus_cols <- grepl("taxa_", names(otus))
rownames(otus) <- otus$public_client_id
sdata <- otus[, !genus_cols]
bowel <- paste0(sdata['bowel'])
table(sdata['bowel'])

   3    1    2    4 
1796   83  769   41 

otus <- as.matrix(otus[, genus_cols])
colnames(otus) <- gsub("taxa_", "", colnames(otus))
tax_matrix <- as.matrix(taxa[, 2:ncol(taxa)])
rownames(tax_matrix) <- taxa[, 1]
as.data.frame(otus)

Now we convert it to fit in a phyloseq object. We need to fulfill the following rules:

1. The OTU table is a matrix with rows being samples and each column being a taxon.
2. The taxonomy table must be a matrix with rows being taxa and columns being ranks.
3. The sample data must be data frame with rows being samples.

To double check if everything works let’s do some validity checks:

stopifnot(all(rownames(otus) %in% rownames(sdata)))
stopifnot(nrow(otus) == nrow(sdata))
print("sample names match")

[1] "sample names match"

stopifnot(all(colnames(otus) %in% rownames(tax_matrix)))
stopifnot(ncol(otus) == nrow(tax_matrix))
stopifnot(!anyDuplicated(tax_matrix))
print("taxa look okay")

[1] "taxa look okay"

stopifnot(!anyDuplicated(sdata))
print("sample data looks okay")

[1] "sample data looks okay"

If that passes we can go ahead and build our phyloseq object.

ps <- phyloseq(
  otu_table(otus, taxa_are_rows = FALSE),
  tax_table(tax_matrix),
  sample_data(sdata)
)

ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 68 taxa and 2689 samples ]
sample_data() Sample Data:       [ 2689 samples by 7 sample variables ]
tax_table()   Taxonomy Table:    [ 68 taxa by 6 taxonomic ranks ]

names(sample_data(ps))

[1] "public_client_id" "sex"              "age"              "BMI_CALC"         "vendor_dashboard" "bowel"            "eGFR"            

bowel <- sample_data(ps)$bowel

#Differential Test
dv_analysis <- differentialTest(formula = ~ bowel + sex + age + BMI_CALC + vendor_dashboard + eGFR,
                                 phi.formula = ~ 1,
                                 formula_null = ~ sex + age + BMI_CALC + vendor_dashboard + eGFR,
                                 phi.formula_null = ~ 1,
                                 data = ps,
                                 test = "LRT", boot = FALSE,
                                 full_output = TRUE,
                                 fdr_cutoff = 0.05)

#See the significant taxa from the computation:
dv_analysis$significant_taxa

 [1] "Bacteria.Actinobacteria.Actinobacteria.Bifidobacteriales.Bifidobacteriaceae.Bifidobacterium"    "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Coriobacteriaceae.Collinsella"         
 [3] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Eggerthellaceae.Adlercreutzia"          "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Bacteroidaceae.Bacteroides"                   
 [5] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Marinifilaceae.Odoribacter"                    "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Rikenellaceae.Alistipes"                      
 [7] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Tannerellaceae.Parabacteroides"                "Bacteria.Firmicutes.Bacilli.Lactobacillales.Streptococcaceae.Streptococcus"                    
 [9] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.Christensenellaceae_R-7_group" "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.nan"                          
[11] "Bacteria.Firmicutes.Clostridia.Clostridiales.Clostridiales_vadinBB60_group.nan"                 "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_AD3011_group"             
[13] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_UCG-001"                   "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.nan"                                  
[15] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Agathobacter"                      "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Blautia"                          
[17] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnoclostridium"                 "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospira"                      
[19] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_NK4A136_group"     "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_UCG-004"          
[21] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Marvinbryantia"                    "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.nan"                              
[23] "Bacteria.Firmicutes.Clostridia.Clostridiales.Peptostreptococcaceae.Romboutsia"                  "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Anaerotruncus"                    
[25] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Butyricicoccus"                    "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Candidatus_Soleaferrea"           
[27] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.DTU089"                            "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Faecalibacterium"                 
[29] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.GCA-900066225"                     "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Intestinimonas"                   
[31] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Negativibacillus"                  "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Oscillibacter"                    
[33] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_5"               "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_9"              
[35] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_NK4A214_group"     "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-002"          
[37] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-004"           "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-005"          
[39] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-013"           "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_1"                   
[41] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_2"                    "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Subdoligranulum"                  
[43] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.UBA1819"                           "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.nan"                              
[45] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.Holdemania"         "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.nan"               
[47] "Bacteria.Verrucomicrobia.Verrucomicrobiae.Verrucomicrobiales.Akkermansiaceae.Akkermansia"      

saveRDS(dv_analysis, "corncob_all3.rds")
#Load the .rds file as "dv_analysis" to continue from already computed dv_analysis file from before.
#This way you don't need to run the computation all over again each time you want to run the code from here below:

library(DataCombine)
taxa<-dv_analysis$all_models
dv_analysis$significant_taxa

 [1] "Bacteria.Actinobacteria.Actinobacteria.Bifidobacteriales.Bifidobacteriaceae.Bifidobacterium"   
 [2] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Coriobacteriaceae.Collinsella"         
 [3] "Bacteria.Actinobacteria.Coriobacteriia.Coriobacteriales.Eggerthellaceae.Adlercreutzia"         
 [4] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Bacteroidaceae.Bacteroides"                   
 [5] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Marinifilaceae.Odoribacter"                   
 [6] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Rikenellaceae.Alistipes"                      
 [7] "Bacteria.Bacteroidetes.Bacteroidia.Bacteroidales.Tannerellaceae.Parabacteroides"               
 [8] "Bacteria.Firmicutes.Bacilli.Lactobacillales.Streptococcaceae.Streptococcus"                    
 [9] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.Christensenellaceae_R-7_group"
[10] "Bacteria.Firmicutes.Clostridia.Clostridiales.Christensenellaceae.nan"                          
[11] "Bacteria.Firmicutes.Clostridia.Clostridiales.Clostridiales_vadinBB60_group.nan"                
[12] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_AD3011_group"             
[13] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.Family_XIII_UCG-001"                  
[14] "Bacteria.Firmicutes.Clostridia.Clostridiales.Family_XIII.nan"                                  
[15] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Agathobacter"                     
[16] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Blautia"                          
[17] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnoclostridium"                
[18] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospira"                      
[19] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_NK4A136_group"    
[20] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Lachnospiraceae_UCG-004"          
[21] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.Marvinbryantia"                   
[22] "Bacteria.Firmicutes.Clostridia.Clostridiales.Lachnospiraceae.nan"                              
[23] "Bacteria.Firmicutes.Clostridia.Clostridiales.Peptostreptococcaceae.Romboutsia"                 
[24] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Anaerotruncus"                    
[25] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Butyricicoccus"                   
[26] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Candidatus_Soleaferrea"           
[27] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.DTU089"                           
[28] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Faecalibacterium"                 
[29] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.GCA-900066225"                    
[30] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Intestinimonas"                   
[31] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Negativibacillus"                 
[32] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Oscillibacter"                    
[33] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_5"              
[34] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminiclostridium_9"              
[35] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_NK4A214_group"    
[36] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-002"          
[37] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-004"          
[38] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-005"          
[39] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcaceae_UCG-013"          
[40] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_1"                   
[41] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Ruminococcus_2"                   
[42] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.Subdoligranulum"                  
[43] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.UBA1819"                          
[44] "Bacteria.Firmicutes.Clostridia.Clostridiales.Ruminococcaceae.nan"                              
[45] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.Holdemania"        
[46] "Bacteria.Firmicutes.Erysipelotrichia.Erysipelotrichales.Erysipelotrichaceae.nan"               
[47] "Bacteria.Verrucomicrobia.Verrucomicrobiae.Verrucomicrobiales.Akkermansiaceae.Akkermansia"      

as.data.frame(dv_analysis$p_fdr)

taxa_p <- DropNA(as.data.frame(dv_analysis$p_fdr), Var = "dv_analysis$p_fdr")

0 rows dropped from the data frame because of missing values.

taxa_p[1]

[1] 0.02120539

taxa[[1]]$coefficients

                                     Estimate  Std. Error     t value     Pr(>|t|)
mu.(Intercept)                   -3.368074680 0.182186294 -18.4869816 6.486060e-72
mu.bowel1                         0.108766498 0.121122817   0.8979852 3.692741e-01
mu.bowel2                         0.136368966 0.046755183   2.9166599 3.567424e-03
mu.bowel4                        -0.212360465 0.178644229  -1.1887340 2.346497e-01
mu.sexW                          -0.212789194 0.048369795  -4.3992163 1.128862e-05
mu.age                           -0.016148200 0.001821779  -8.8639749 1.378229e-18
mu.BMI_CALC                      -0.003154924 0.003370311  -0.9360928 3.493098e-01
mu.vendor_dashboardSecond Genome  0.127694322 0.060993339   2.0935782 3.639137e-02
mu.eGFR                           0.001582766 0.001243639   1.2726892 2.032388e-01
phi.(Intercept)                  -2.808912916 0.043294681 -64.8789385 0.000000e+00

sample_data(arivale_phylo) 

Sample Data:        [3653 samples by 26 sample variables]:

#Tom's code:
#Rarefy genotek dataset to even depth
#rarefy even depth is a command in phyloseq
#rarefy for alpha diversity, but not for corncob
rarefied_genotek=rarefy_even_depth(arivale_phylo, sample.size = min(sample_sums(arivale_phylo)),
  rngseed = 111, replace = FALSE, trimOTUs = TRUE, verbose = TRUE)

`set.seed(111)` was used to initialize repeatable random subsampling.
Please record this for your records so others can reproduce.
Try `set.seed(111); .Random.seed` for the full vector
...
25984OTUs were removed because they are no longer 
present in any sample after random subsampling

...

richness <- estimate_richness(rarefied_genotek, measures=c("Shannon","Observed"))
#Refer to Christian (phytree?) on PD whole tree (phyloseq additional PD whole tree package, check it) and add it to richness column on df.

richness$public_client_id <-sample_data(rarefied_genotek)$public_client_id
#Refer to Christian (phytree?) on PD whole tree (phyloseq additional PD whole tree package, check it) and add it to richness column on df.

library(btools)
richness$Pielou <- richness$Shannon/richness$Observed

richness

hist(richness$Pielou,
     main = "Pielou's Evenness (Shannon/Observed ASV)",
     xlab = "Evenness",
     xlim = c(0,0.05),
     breaks = 200)

hist(richness$Shannon,
     main = "Shannon Diversity",
     xlab = "Diversity Index",
     xlim = c(0,6),
     breaks = 20)

hist(richness$Observed,
     main = "Observed ASVs",
     xlab = "ASVs",
     xlim = c(0,1000),
     breaks = 30)

library(tidyverse)
library(ggpmisc)
library(quantreg)
library(broom)
library(ggbeeswarm)
library(ggsignif)
library(gginnards)
library(ggpubr)
library(ggpmisc)
library(broom.mixed)
library(scales)
#library(nlme)

sam.n <- function(x){
  return(c(y = mean(x), label = length(x)))
}

combinations <- list(c("Constipation","Low Normal"),
                                 c("Constipation","High Normal"),
                                 c("Constipation","Diarrhea"),
                                 c("Low Normal","High Normal"),
                                 c("Low Normal","Diarrhea"),
                                 c("High Normal","Diarrhea"))
dat$bowel <- factor(dat$bowel, levels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))

library(compositions)
df <- read.csv("gut_full.csv", check.names=F)
colnames(df) = gsub("nan", "Unclassified", colnames(df))
comparisons = list(c("Constipation","High Normal"),c("Low Normal","High Normal"),c("Diarrhea","High Normal"))
df$bowel <- factor(df$bowel, levels = c(1,2,3,4), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))
df_otus <- dplyr::select(df, -c("public_client_id","bowel","vendor_dashboard","sex","age","BMI_CALC","eGFR"))

df_otus <- as.data.frame(clr(as.matrix(df_otus)))
df_select <- dplyr::select(df, c(1:7))
df_otus <- cbind(df_select,df_otus)
df <- df_otus
dat <- df_otus

library(ggrepel)
sam.n <- function(x){return(c(y = mean(x), label = length(x)))}

# ggplot2 plots
library(gplots)
library(ggrepel)
library(ggbreak)
library(sommer)

#BMF df
const <- c()
const_p <- c()
low <- c()
low_p <- c()
diarrhea <- c()
diarrhea_p <- c()
taxa_names <- c()
family_names <- c()
genus_names <- c()
likelihood <- c()
adj_p <- c()
names(dv_analysis$p) = gsub("nan", "Unclassified", names(dv_analysis$p))

#dv_analysis$all_models[order(dv_analysis$p[order(names(dv_analysis$p))])])
for (i in 1:length(dv_analysis$all_models[names(dv_analysis$p)])) {
  const_p[i] <- dv_analysis$all_models[[i]]$coefficients[2,4]
  const[i] <- dv_analysis$all_models[[i]]$coefficients[2,1]
  low_p[i] <- dv_analysis$all_models[[i]]$coefficients[3,4]
  low[i] <- dv_analysis$all_models[[i]]$coefficients[3,1]
  diarrhea_p[i] <- dv_analysis$all_models[[i]]$coefficients[4,4]
  diarrhea[i] <- dv_analysis$all_models[[i]]$coefficients[4,1]
  taxa_names[i] <- names(dv_analysis$p)[i]
  family_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][5]
  genus_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][6]
  likelihood[i] <- dv_analysis$all_models[[i]]$logL
  adj_p[i] <- dv_analysis$p_fdr[[i]]
}

genus_names[i] <- strsplit(taxa_names[i],'.',fixed=TRUE)[[1]][6]

p_df <- bind_cols(taxa_names,family_names,genus_names, likelihood, adj_p, const,const_p,low,low_p,diarrhea,diarrhea_p)

New names:
• `` -> `...1`
• `` -> `...2`
• `` -> `...3`
• `` -> `...4`
• `` -> `...5`
• `` -> `...6`
• `` -> `...7`
• `` -> `...8`
• `` -> `...9`
• `` -> `...10`
• `` -> `...11`

names(p_df) <- c("Genera","Family","Genus","LogL","Adj.P","Const.Beta","Const.P.Val","Low.Beta","Low.P.Val","Diarrhea.Beta","Diarrhea.P.Val")
p_df$Combined <- paste(p_df$Family,p_df$Genus)
ab <- colSums(df_otus[names(df_otus) %in% paste("taxa_",names(dv_analysis$p),sep="")])/colSums(!!df_otus[names(df_otus) %in% paste("taxa_",names(dv_analysis$p),sep="")])
p_df$Mean.Abundance <- ab[paste("taxa_",p_df$Genera,sep="")]

p_df$Const.Adj.P.Val <- p.adjust(const_p, method = "fdr", n = length(const_p))
p_df$Low.Adj.P.Val <- p.adjust(low_p, method = "fdr", n = length(low_p))
p_df$Diarrhea.Adj.P.Val <- p.adjust(diarrhea_p, method = "fdr", n = length(diarrhea_p))

p_df <- p_df[order(-p_df$Mean.Abundance),]


set <- subset(p_df[which(p_df$Adj.P < 0.05),], Genus != 'Unclassified')
set <- subset(set[order(-set$Mean.Abundance),], Genus != 'Unclassified')


abundant_list <- list()
abundant_genus <- list()
abundant <- paste("taxa_",set[order(-set$Mean.Abundance),]$Genera,sep="")
for (i in 1:68) {
  abundant_list[i] <- paste(strsplit(abundant[i],'.',fixed=TRUE)[[1]][5],strsplit(abundant[i],'.',fixed=TRUE)[[1]][6])
  abundant_genus[i] <- paste(strsplit(abundant[i],'.',fixed=TRUE)[[1]][6])
}

#top 10 most abundant non-unclassified hits
abundant_trunc <- set[match(abundant_list[!is.na(set$Combined[match(abundant_list,set$Combined)])],set$Combined),][which(set$Adj.P<0.05),]$Combined[1:10]
`%!in%` <- Negate(`%in%`)
t <- set[set$Combined %!in% abundant_trunc,]
#t[!is.na(t[order(t$Adj.P),]$Mean.Abundance),][order(t$Adj.P),]

span <- rbind(set[order(set[-order(set[!set$Combined %in% abundant_trunc,]$Adj.P),][which(set$Adj.P<0.05),]$Mean.Abundance %!in% abundant_trunc),][1:10,],
t[match('Akkermansiaceae Akkermansia',t$Combined),],
t[order(t$Const.Adj.P.Val),][1:9,])
span$Letters <- LETTERS[1:20]

p_df$Letters <- c(NA)
test <- rbind(span,p_df[p_df$Combined %!in% span$Combined,])
test <- test[!duplicated(test$Genera),]
#test[order(test[!duplicated(test),]$Letters),]
p_df <- test
p_df <- p_df[order(p_df$Adj.P),]
p_df <- p_df[order(-p_df$Mean.Abundance),]


model = y ~ x
df_num <- merge(richness, df_otus, by="public_client_id")
df_num <- within(df_num, bowel <- relevel(bowel, ref = 'High Normal'))
d1_plot <- lm(Pielou ~ bowel, data = df_num)
pval1 = summary(d1_plot)$coefficients[2,4]
r1 = summary(d1_plot)$adj.r.squared
pval1

[1] 0.001840045

r1

[1] 0.01204487

d1 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Pielou, group = factor(bowel))) +
  scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r1),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval1),3)))))+
  geom_beeswarm(aes(color = bowel), cex = 0.5) +
  geom_boxplot(alpha=0) +
  ylim(0.005,0.05)+
  geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
  ggtitle(label = "C)\nPielou's Evenness vs BMF") +
  xlab("Bowel Movement Frequency (BMF)") +
  ylab("Pielou's Evenness") +
  scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
  guides(colour = guide_legend(override.aes = list(size=7))) +
  theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d1

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (position_beeswarm).

model = y ~ x
d2_plot <- lm(Shannon ~ bowel, data = df_num)
pval2 = summary(d2_plot)$coefficients[2,4]
r2 = summary(d2_plot)$adj.r.squared
pval2

[1] 0.00964909

r2

[1] 0.01384

d2 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Shannon, group = bowel)) +
  scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r2),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval2),3)))))+
  geom_beeswarm(aes(color = factor(bowel)), cex = 0.5) +
  geom_boxplot(alpha=0) +
  geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
  ggtitle(label = "B)\nShannon Diversity vs BMF") +
  xlab("Bowel Movement Frequency (BMF)") +
  ylab("Shannon Diversity") +
  scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
  guides(colour = guide_legend(override.aes = list(size=7))) +
  theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d2

model = y ~ x
d3_plot <- lm(Observed ~ bowel, data = df_num)
pval3 = summary(d3_plot)$coefficients[2,4]
r3 = summary(d3_plot)$adj.r.squared
pval3

[1] 5.448268e-05

r3

[1] 0.0208449

d3 <- ggplot(data = df_num, aes(x = factor(bowel, level = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea")), y = Observed, group = bowel))+
  scale_x_discrete(bquote(atop(~italic("adj R")^2~" = "~.(formatC(as.numeric(r3),3)),~italic("P value")~" = "~.(formatC(as.numeric(pval3),3)))))+
  geom_beeswarm(aes(color = factor(bowel)), cex = 0.5) +
  geom_boxplot(alpha=0) +
  geom_smooth(method = "lm", formula = y ~ x, aes(group = 1)) +
  ggtitle(label = "A)\nObserved ASVs vs BMF") +
  xlab("Bowel Movement Frequency (BMF)") +
  ylab("Observed ASVs") +
  scale_colour_discrete(name="BMF Category", limits = c("Constipation", "Low Normal", "High Normal", "Diarrhea"), labels=c("Constipation", "Low Normal", "High Normal", "Diarrhea")) +
  guides(colour = guide_legend(override.aes = list(size=7))) +
  theme(plot.title = element_text(size=10), legend.text = element_text(size = 10), legend.title = element_text(size = 10), axis.text.x = element_blank(), axis.title.y = element_text(size = 10))
d3

PSO <- ggarrange(d3, d2, d1, legend = "top", align = "hv", common.legend = TRUE, widths = c(7,7,7), heights = c(20,20,20), nrow = 1, ncol = 3)+ theme(plot.margin = margin(0,0,0,0))

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (position_beeswarm).

PSO
ggsave(
  "PSOvBMF.png",
  plot = PSO,
  device = NULL,
  path = NULL,
  scale = 1,
  width = NA,
  height = NA,
  units = c("in", "cm", "mm", "px"),
  dpi = 300,
  limitsize = TRUE,
  bg = NULL
)

Saving 7.29 x 4.51 in image

# Significant Genera Plot (plot 0)
plot0s <- ggplot(p_df, mapping = aes(x=Const.Beta,y = -log10(Const.P.Val), color = ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"),group = factor(ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"))))+
  geom_vline(xintercept = 0)+
  geom_point(size=1, aes(x = Const.Beta, y = -log10(Const.P.Val), group = factor(ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Letters,""))),position = position_dodge(width = 1)) +
  geom_label_repel(aes(x = Const.Beta, y = -log10(Const.P.Val), group = factor(ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Letters,""))), label.size = 0.1, label.padding = 0.1, box.padding = 0.5, min.segment.length = 0, point.size = NA, position = position_dodge(width = 1), size=2,label = factor(ifelse(p_df$Const.Adj.P.Val < 0.05,p_df$Letters,"")),color = "black", check_overlap = FALSE, max.overlaps = 127)+
 #ggtitle("Significant Genera") +
  xlab(bquote(beta~" Coefficient")) + 
  ylab(bquote("-log"[10]~"(P value)")) + 
  #scale_x_break(c(-0.6, -18)) +
  scale_x_continuous(name = bquote(atop(beta["BMF"]~" Coefficient",italic("Constipation"))), 
                     guide = guide_axis(n.dodge = 2), limits = c(-1.5,1.5)) +
  theme(text = 
          element_text(size = 14), 
        plot.title = element_text(vjust = 0.5), 
        #plot.subtitle = element_text(size=8, hjust = 0.5), 
        legend.title = element_blank(), 
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 8),
              legend.text  = element_text(size = 3),
              legend.position = "top",
              legend.key.size = unit(0.1, "lines"))+
          guides(shape = guide_legend(override.aes = list(size = 1.5)),
               color = guide_legend(override.aes = list(size = 1.5), nrow = 25),
               fill=guide_legend(title=NULL),
               aspect.ratio = 0.95)

Warning: Ignoring unknown parameters: check_overlap

plot0s

Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
Warning: Removed 13 rows containing missing values (geom_label_repel).

# Significant Genera Plot (plot 0)
plot0m <- ggplot(p_df, mapping = aes(x=Low.Beta,y = -log10(Low.P.Val), color = ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"),group = factor(ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"))))+
  geom_vline(xintercept = 0)+
  geom_point(size=1, aes(x = Low.Beta, y = -log10(Low.P.Val), group = factor(ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Letters,""))),position = position_dodge(width = 1)) +
  geom_label_repel(aes(x = Low.Beta, y = -log10(Low.P.Val), group = factor(ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Letters,""))), label.size = 0.1, label.padding = 0.1, box.padding = 0.5, min.segment.length = 0, point.size = NA, position = position_dodge(width = 1), size=2,label = factor(ifelse(p_df$Low.Adj.P.Val < 0.05,p_df$Letters,"")),color = "black", check_overlap = FALSE, max.overlaps = 127)+
 #ggtitle("Significant Genera") +
  xlab(bquote(beta~" Coefficient")) + 
  ylab(bquote("-log"[10]~"(P value)")) +
  #scale_x_break(c(-0.6, -18)) +
  scale_x_continuous(name = bquote(atop(beta["BMF"]~" Coefficient",italic("Low Normal"))), 
                     guide = guide_axis(n.dodge = 2)) +
  theme(text = 
          element_text(size = 14), 
        plot.title = element_text(vjust = 0.5), 
        #plot.subtitle = element_text(size=8, hjust = 0.5), 
        legend.title = element_blank(), 
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 8),
              legend.text  = element_text(size = 3),
              legend.position = "top",
              legend.key.size = unit(0.1, "lines"))+
          guides(shape = guide_legend(override.aes = list(size = 1.5)),
               color = guide_legend(override.aes = list(size = 1.5), nrow = 25),
               fill=guide_legend(title=NULL),
               aspect.ratio = 0.95)

Warning: Ignoring unknown parameters: check_overlap

plot0m

Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
Warning: Removed 27 rows containing missing values (geom_label_repel).

# Significant Genera Plot (plot 0)
plot0d <- ggplot(p_df, mapping = aes(x=Diarrhea.Beta,y = -log10(Diarrhea.P.Val), color = ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"),group = factor(ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Genera,"- FDR Adj. P Value"))))+
  geom_vline(xintercept = 0)+
  geom_point(size=1, aes(x = Diarrhea.Beta, y = -log10(Diarrhea.P.Val), group = factor(ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Letters,""))),position = position_dodge(width = 1)) +
  geom_label_repel(aes(x = Diarrhea.Beta, y = -log10(Diarrhea.P.Val), group = factor(ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Letters,""))), label.size = 0.1, label.padding = 0.1, box.padding = 0.5, min.segment.length = 0, point.size = NA, position = position_dodge(width = 1), size=2,label = factor(ifelse(p_df$Diarrhea.Adj.P.Val < 0.05,p_df$Letters,"")),color = "black", check_overlap = FALSE, max.overlaps = 127)+
 #ggtitle("Significant Genera") +
  xlab(bquote(beta~" Coefficient")) + 
  ylab(bquote("-log"[10]~"(P value)")) +
  #scale_x_break(c(-0.6, -18)) +
  scale_x_continuous(name = bquote(atop(beta["BMF"]~" Coefficient",italic("Diarrhea"))), 
                     guide = guide_axis(n.dodge = 2)) +
  theme(text = 
          element_text(size = 14), 
        plot.title = element_text(vjust = 0.5), 
        #plot.subtitle = element_text(size=8, hjust = 0.5), 
        legend.title = element_blank(), 
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 8),
              legend.text  = element_text(size = 3),
              legend.position = "top",
              legend.key.size = unit(0.1, "lines"))+
          guides(shape = guide_legend(override.aes = list(size = 1.5)),
               color = guide_legend(override.aes = list(size = 1.5), nrow = 25),
               fill=guide_legend(title=NULL),
               aspect.ratio = 0.95)

Warning: Ignoring unknown parameters: check_overlap

plot0d

Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals

library(ggpubr)
figure1 <- ggarrange(plot0s,plot0m,plot0d, common.legend = FALSE, nrows = 1, widths=1,heights=3,ncols = 3, align = "hv", legend = "top") + theme(plot.margin = margin(0,0,0,0))

Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Const.Adj.P.Val` is discouraged. Use `Const.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
Warning: Removed 13 rows containing missing values (geom_label_repel).
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Low.Adj.P.Val` is discouraged. Use `Low.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals
Warning: Removed 27 rows containing missing values (geom_label_repel).
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Genera` is discouraged. Use `Genera` instead.
Warning: Use of `p_df$Diarrhea.Adj.P.Val` is discouraged. Use `Diarrhea.Adj.P.Val` instead.
Warning: Use of `p_df$Letters` is discouraged. Use `Letters` instead.
Warning: position_dodge requires non-overlapping x intervals
Warning: position_dodge requires non-overlapping x intervals

ggsave(
  "BMFvsGenera1.png",
  plot = figure1,
  device = NULL,
  path = NULL,
  scale = 1.75,
  width = NA,
  height = NA,
  units = c("in", "cm", "mm", "px"),
  dpi = 300,
  limitsize = TRUE,
  bg = NULL
)

Saving 12.8 x 7.88 in image

figure1

Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  : 
  Viewport has zero dimension(s)

sig = function(x){
  if(x < 0.001){"***"} 
  else if(x < 0.01){"**"}
  else if(x < 0.05){"*"}
  else{NA}}

df <- read.csv("gut_full.csv", check.names=F)
colnames(df) = gsub("nan", "Unclassified", colnames(df))
comparisons = list(c("Constipation","High Normal"),c("Low Normal","High Normal"),c("Diarrhea","High Normal"))
df$bowel <- factor(df$bowel, levels = c(1,2,3,4), labels = c("Constipation", "Low Normal", "High Normal", "Diarrhea"))
df_otus <- dplyr::select(df, -c("public_client_id","bowel","vendor_dashboard","sex","age","BMI_CALC","eGFR"))

df_otus <- as.data.frame(clr(as.matrix(df_otus)))
df_select <- dplyr::select(df, c(1:7))
df_otus <- cbind(df_select,df_otus)
df <- df_otus
dat <- df_otus

#Gut Genera
df$p.value <- NA
df$pval <- NULL
df$p.val <- NULL


counter <<- 1

test = function(x,y,z,j) {
  work <- NULL
  results <- NULL
  x = NULL
  y = NULL
  if (counter == 1) {z = comparisons[[1]][1]}
  else if (counter == 2) {z = comparisons[[2]][1]}
  else if (counter > 2) {z = comparisons[[3]][1]}
  temp <- data.frame("name" = NA, "p.value" = NA, "bowel" = NA)
  for (i in 1:20) {
    temp$name <- paste("genus_",span$Genera[[i]],sep="")
    temp$p.value <- span$Const.Adj.P.Val[[i]] 
    temp$bowel <- "Constipation"
    work <- rbind(work,temp)
    temp$p.value <- span$Low.Adj.P.Val[[i]] 
    temp$bowel <- "Low Normal"
    work <- rbind(work,temp)
    temp$p.value <- span$Diarrhea.Adj.P.Val[[i]] 
    temp$bowel <- "Diarrhea"
    work <- rbind(work,temp)
  }
  if (z[[1]][1] == comparisons[[1]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1])}
  else if (z[[1]][1] == comparisons[[2]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1+1]) }
  else if (z[[1]][1] == comparisons[[3]][1]) {results <- list(p.value = work$p.value[3*(j-1)+1+2])}
  counter<<-counter+1
  if (counter > 3) {counter<<-1}
  return(results)
}

myplots <- list()  # new empty list
df_test <- df
df_test[df_test==0] <- NA
for (genera in 1:20) {
  y_name <- paste("taxa_",span$Genera[genera],sep="")
  myplots[[genera]] <- local({
    plotlim_lower = min(
df_test[!is.na(df_test[y_name]),][y_name])
    plotlim_upper = max(
df_test[!is.na(df_test[y_name]),][y_name])
    plotlim_bar = plotlim_lower - 4
    plotlim_margin = 6
    genera <- genera
    plt <- ggplot(data = df_test, aes(x = bowel, y = .data[[y_name]], group = bowel)) +
    scale_x_discrete(guide = guide_axis(n.dodge = 2))+
    geom_jitter(aes(color = bowel),  size = 0.1, cex = 0.05) +
    geom_boxplot(data = df_test, alpha=0.0,outlier.shape = NA) +
    theme(text = element_text(size = 9)) +
    ggtitle(label = paste(span$Family[genera],"\n",span$Genus[genera],sep="")) +
    geom_signif(data = df_test, comparisons = comparisons, map_signif_level = sig, test = 'test', test.args = list(z = comparisons, j = genera), y_position = plotlim_bar,step_increase = 0.15,  size = 0.5, 
                textsize = 1.5,
                tip_length = c(0,0)) +
    coord_cartesian(ylim=c(plotlim_lower,plotlim_upper),clip="off")+
    labs(color = "BMF Category", y = ifelse((genera == 1 | genera  == 6 | genera  == 11 | genera == 16),"CLR Abundance","")) +
    guides(colour = guide_legend(override.aes = list(size=7), title.position = 'left', nrow = 1, ncol = 4)) +
    theme(plot.margin = unit(c(0,0,plotlim_margin,0), "cm"), 
          plot.title = element_text(size=5.75), 
          legend.title = element_text(size=10), 
          plot.subtitle = element_text(size=10), 
          legend.text = element_text(size=7),
          axis.text.x = element_blank(), 
          axis.text.y = element_text(size=7), 
          axis.title.y = element_text(size=7),
          axis.title.x = element_blank(),
          aspect.ratio = 0.95)
  })
}

Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size
Warning: Duplicated aesthetics after name standardisation: size

library(gplots)
library(ggrepel)
library(ggbreak)

#Truncate the top 20 including verukkamansia 
title <- list()
for(i in seq(1:20)) { title[i] <- myplots[[c(i)]]$labels$title }
plots1 <- myplots[title %in% paste(
  span$Family,
  "\n",
  span$Genus,sep="")][1:20]
plots_final <- plots1[!sapply(plots1,is.null)]

final1 <- ggarrange(plotlist = plots_final[1:5], labels = LETTERS[1:5], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))

Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 56 rows containing non-finite values (stat_boxplot).
Warning: Removed 56 rows containing non-finite values (stat_signif).
Warning: Removed 56 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing non-finite values (stat_signif).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 85 rows containing non-finite values (stat_boxplot).
Warning: Removed 85 rows containing non-finite values (stat_signif).
Warning: Removed 85 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).

final1

final2 <- ggarrange(plotlist = plots_final[6:10], labels = LETTERS[6:10], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))

Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 98 rows containing non-finite values (stat_boxplot).
Warning: Removed 98 rows containing non-finite values (stat_signif).
Warning: Removed 98 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 607 rows containing non-finite values (stat_boxplot).
Warning: Removed 607 rows containing non-finite values (stat_signif).
Warning: Removed 607 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 4 rows containing non-finite values (stat_boxplot).
Warning: Removed 4 rows containing non-finite values (stat_signif).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 7 rows containing missing values (geom_signif).

final2

final3 <- ggarrange(plotlist = plots_final[11:15], labels = LETTERS[11:15], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))

Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 271 rows containing non-finite values (stat_boxplot).
Warning: Removed 271 rows containing non-finite values (stat_signif).
Warning: Removed 271 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 591 rows containing non-finite values (stat_boxplot).
Warning: Removed 591 rows containing non-finite values (stat_signif).
Warning: Removed 591 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 147 rows containing non-finite values (stat_boxplot).
Warning: Removed 147 rows containing non-finite values (stat_signif).
Warning: Removed 147 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 240 rows containing non-finite values (stat_boxplot).
Warning: Removed 240 rows containing non-finite values (stat_signif).
Warning: Removed 240 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).

final3

final4 <- ggarrange(plotlist = plots_final[16:20], labels = LETTERS[16:20], legend = "top", align = "hv", font.label = list(size = 8), common.legend = TRUE, nrow = 1, ncol = 5) + theme(plot.margin = unit(c(0,0,0,0), "cm"))

Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 40 rows containing non-finite values (stat_boxplot).
Warning: Removed 40 rows containing non-finite values (stat_signif).
Warning: Removed 40 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 331 rows containing non-finite values (stat_boxplot).
Warning: Removed 331 rows containing non-finite values (stat_signif).
Warning: Removed 331 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 680 rows containing non-finite values (stat_boxplot).
Warning: Removed 680 rows containing non-finite values (stat_signif).
Warning: Removed 680 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 20 rows containing non-finite values (stat_boxplot).
Warning: Removed 20 rows containing non-finite values (stat_signif).
Warning: Removed 20 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).

final4

counter<<-1
figure_main <- ggarrange(plotlist = plots_final[1:20], labels = LETTERS[1:20], legend = "top", align = "hv", common.legend = TRUE, nrow = 4, ncol = 5)

Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 3 rows containing non-finite values (stat_boxplot).
Warning: Removed 3 rows containing non-finite values (stat_signif).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 56 rows containing non-finite values (stat_boxplot).
Warning: Removed 56 rows containing non-finite values (stat_signif).
Warning: Removed 56 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Warning: Removed 2 rows containing non-finite values (stat_signif).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 85 rows containing non-finite values (stat_boxplot).
Warning: Removed 85 rows containing non-finite values (stat_signif).
Warning: Removed 85 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 113 rows containing non-finite values (stat_boxplot).
Warning: Removed 113 rows containing non-finite values (stat_signif).
Warning: Removed 113 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 98 rows containing non-finite values (stat_boxplot).
Warning: Removed 98 rows containing non-finite values (stat_signif).
Warning: Removed 98 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 607 rows containing non-finite values (stat_boxplot).
Warning: Removed 607 rows containing non-finite values (stat_signif).
Warning: Removed 607 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 171 rows containing non-finite values (stat_boxplot).
Warning: Removed 171 rows containing non-finite values (stat_signif).
Warning: Removed 171 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 4 rows containing non-finite values (stat_boxplot).
Warning: Removed 4 rows containing non-finite values (stat_signif).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 7 rows containing missing values (geom_signif).
Warning: Removed 762 rows containing non-finite values (stat_boxplot).
Warning: Removed 762 rows containing non-finite values (stat_signif).
Warning: Removed 762 rows containing missing values (geom_point).
Warning: Removed 8 rows containing missing values (geom_signif).
Warning: Removed 271 rows containing non-finite values (stat_boxplot).
Warning: Removed 271 rows containing non-finite values (stat_signif).
Warning: Removed 271 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 591 rows containing non-finite values (stat_boxplot).
Warning: Removed 591 rows containing non-finite values (stat_signif).
Warning: Removed 591 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 147 rows containing non-finite values (stat_boxplot).
Warning: Removed 147 rows containing non-finite values (stat_signif).
Warning: Removed 147 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 240 rows containing non-finite values (stat_boxplot).
Warning: Removed 240 rows containing non-finite values (stat_signif).
Warning: Removed 240 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 592 rows containing non-finite values (stat_boxplot).
Warning: Removed 592 rows containing non-finite values (stat_signif).
Warning: Removed 592 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 40 rows containing non-finite values (stat_boxplot).
Warning: Removed 40 rows containing non-finite values (stat_signif).
Warning: Removed 40 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 331 rows containing non-finite values (stat_boxplot).
Warning: Removed 331 rows containing non-finite values (stat_signif).
Warning: Removed 331 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 680 rows containing non-finite values (stat_boxplot).
Warning: Removed 680 rows containing non-finite values (stat_signif).
Warning: Removed 680 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).
Warning: Removed 20 rows containing non-finite values (stat_boxplot).
Warning: Removed 20 rows containing non-finite values (stat_signif).
Warning: Removed 20 rows containing missing values (geom_point).
Warning: Removed 5 rows containing missing values (geom_signif).

ggsave(
  "BMFvsGenera2.png",
  plot = final1,
  device = NULL,
  path = NULL,
  scale = 1.5,
  width = NA,
  height = NA,
  units = c("in", "cm", "mm", "px"),
  dpi = 300,
  limitsize = TRUE,
  bg = NULL
)

Saving 10.9 x 6.76 in image

ggsave(
  "BMFvsGenera3.png",
  plot = final2,
  device = NULL,
  path = NULL,
  scale = 1.5,
  width = NA,
  height = NA,
  units = c("in", "cm", "mm", "px"),
  dpi = 300,
  limitsize = TRUE,
  bg = NULL
)

Saving 10.9 x 6.76 in image

ggsave(
  "BMFvsGenera4.png",
  plot = final3,
  device = NULL,
  path = NULL,
  scale = 1.5,
  width = NA,
  height = NA,
  units = c("in", "cm", "mm", "px"),
  dpi = 300,
  limitsize = TRUE,
  bg = NULL
)

Saving 10.9 x 6.76 in image

ggsave(
  "BMFvsGenera5.png",
  plot = final4,
  device = NULL,
  path = NULL,
  scale = 1.5,
  width = NA,
  height = NA,
  units = c("in", "cm", "mm", "px"),
  dpi = 300,
  limitsize = TRUE,
  bg = NULL
)

Saving 10.9 x 6.76 in image

ggsave(
  "BMFvsGeneraMain.png",
  plot = final_main,
  device = NULL,
  path = NULL,
  scale = 1.5,
  width = NA,
  height = NA,
  units = c("in", "cm", "mm", "px"),
  dpi = 300,
  limitsize = TRUE,
  bg = NULL
)

Saving 10.9 x 6.76 in image

ggsave(
  "BMFvsGeneraMain.png",
  plot = final_main,
  device = NULL,
  path = NULL,
  scale = 0.5,
  width = NA,
  height = NA,
  units = c("in", "cm", "mm", "px"),
  dpi = 300,
  limitsize = TRUE,
  bg = NULL
)

Saving 7 x 7 in image